package scu.maqiang.heat;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.*;

public class ConvectDiffusion2DT31 {
    /**
     * @param args
     */
    public static void main(String[] args) {
        long startTime = System.nanoTime();
        Mesh2T3 Th = new Mesh2T3().square2D(100, 100, T3Type.Right);
        Th.scale(2, 2);
        Th.translate(-1, -1);
        //Th.scale(1, Math.PI * 3 / 4);
        //Th.transform(ff);
        //Th.writeToTecplot("mesh.dat");

        FES2T31 fespace = new FES2T31(Th);
        SRMatrix A = new SRMatrix(fespace.GetNdof());
        fespace.assembleStiff(new double[]{0.1}, BVPType.COMMON, A);
        //System.out.println("Number of A " + A.size());
        //System.out.println("A = " + MMatrix.toString(A.full()));
//
//          MQSparseMatrix L, U;
//          A.LUFacterization(out L, out U);

        //Console.WriteLine(A.getNumber() + "  " + L.getNumber() + "  " + U.getNumber());
        //Console.WriteLine("L = " + L);
        //Console.WriteLine("U = " + U);

        //Console.WriteLine("L*U = " + L.full()*U.full());
//
        fespace.assembleConvection(new double[]{1.0, 1.0}, BVPType.COMMON, A);
        A.sort();
        //System.out.println(A);
//          //MQMatrix fullB = B.full();
//          //Console.WriteLine("B = " + fullB);
//
        //System.out.println("B = " + MMatrix.toString(B.full()));
        //SRMatrix AA = A.add(B);
        //ConstFunction source = new ConstFunction(10);
        //System.out.println("AA = " + MMatrix.toString(AA.full()));
        double[] source = new double[fespace.GetNdof()];
        fespace.assembleSource(new double[]{1.0}, BVPType.COMMON, source);
        //System.out.println(MVO.toString(source));

        //double[] xx = new double[fespace.GetNdof()];
        // Arrays.fill(xx, 1.0);
        //System.out.println(MVO.toString(A.mul(xx)));
//
//          //ConstFunction g = new ConstFunction(10);
//          //MQVector flux = fespace.assembleHeatFlux(g, 1, 1);
//
        //fespace.applyBC_MBN(A, Direct.X, 1, 2, 3, 4);
        //fespace.applyBC_MBN(source, Direct.X, 0.0, 1, 2, 3, 4);
        fespace.applyBCZero(A, source, Direct.All, BCType.RCE, 1, 2, 3, 4);
//        NewIterSSolver solver = new NewIterSSolver(A);
//        solver.setEps(1.0e-10);
//
//        //IterSSolver solver = new IterSSolver();
//        //A.sort();
//        //System.out.println(A);
//        //System.out.println(MVO.toString(source));
//
//        double[] x = new double[fespace.GetNdof()];
//        solver.GMRES_SSOR(source, x, 1.5, 10, 1);
//        System.out.println("min: " + MVO.min(x) + "  max: " + MVO.max(x));
//
//        Arrays.fill(x, 0.0);
//        solver.PCGSSSOR(source, x, 1.5, 1);
//        System.out.println("min: " + MVO.min(x) + "  max: " + MVO.max(x));
//
//        Arrays.fill(x, 0.0);
//        solver.CGNE(source, x, 1);
//        System.out.println("min: " + MVO.min(x) + "  max: " + MVO.max(x));

        DirectSSolver directSolver = new DirectSSolver(A);
        double[] x = directSolver.LUSolve(source, 1, 1e-6);
//        System.out.println("min: " + MVO.min(x) + "  max: " + MVO.max(x));

//         double[] x1 = new double[fespace.getNDof()];
//          solver.CGNR(AA, RHS, x1);
//          System.out.println("min: " + MVector.min(x1) + "  max: " + MVector.max(x1));
//          
//          System.out.println("Error = " + MVector.distance(x0, x1)/MVector.L2Norm(x1));
        Th.toTecplot("solution2.dat", x);
        //System.out.println(MVO.toString(x));
        long endTime = System.nanoTime();
        System.out.println("Elapse Time: " + (endTime - startTime) / 1e9 + " s");
    }

}
